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SUMMARY 


A  Muskingum-Cunge  channel  flow  routing  scheme  is  modified  for  application  to 
large  drainage  networks  with  compound  cross  sections  and  for  continuous  long-term 
simulation.  The  modifications  consist  of  a  decoupling  and  separate  routing  of  main  and 
overbank  channel  flow,  an  introduction  of  a  variable  time  step  to  increase  model  efficiency 
during  periods  of  steady  flow,  and  an  internal  determination  of  the  numerical  increment. 
The  resulting  hydrologic  model  is  verified  by  comparing  its  flow  routing  results  with  those 
of  hydraulic  benchmark  models  solving  the  full  unsteady  flow  equations.  Test  conditions 
consist  of  hypothetical  flood  hydrographs,  long  prismatic  channels  with  simple  and 
compound  sections,  and  a  third  order  drainage  network.  For  the  tested  conditions,  the 
model  produces  hydrograph  peaks,  times  to  peak  and  shapes  that  compare  well  with 
those  of  the  hydraulic  benchmarks.  Hydrograph  distortions  due  to  overbank  flood  plain 
storage  and  multiple  peaks  from  complex  drainage  networks  are  also  well  reproduced. 
The  execution  time  of  the  model  is  generally  one  order  of  magnitude  faster  than  that  of 
the  hydraulic  benchmark  models. 


Copy  1 

■'‘s^ccrto^ 

1  Aoasaslon  for  /'  ] 

FTIS  ORAJtl 

nr 

DTIC  TAB 

□ 

Unanno’Jinced 

□ 

Justification- 

___ - 

niRtributlOO/ 

Aval  lab  ilU7 

Cedes 

Dl8t 


lAvail  and/or 
Special 


3 


INTRODUCTION 

In  drainage  networks,  channel  runoff  is  continuously  transformed  as  it  travels 
downstream.  Within  individual  channel  reaches,  hydrograph  characteristics  change  as 
a  result  of  flow  hydraulics,  channel  storage,  subsurface  contributions  or  transmission 
losses,  and  lateral  inflow.  In  the  presence  of  active  flood  plains,  overbank  storage 
produces  additional  attenuation.  And,  at  channel  junctions  discreet  changes  in  runoff 
occur  as  a  result  of  the  merging  of  flows  from  upstream  areas.  These  flow 
transformations  occur  simultaneously  throughout  the  drainage  network  and  reshape  the 
channel-flow  hydrographs  as  they  travel  downstream.  In  addition  to  these  in-channel 
transformations,  the  spatial  distribution  of  the  source  areas  and  the  timing  of  their 
respective  runoff  stongly  influence  the  temporal  characteristics  of  the  watershed  response. 
To  quantify  these  effects  in  large  watersheds  and  complex  drainage  networks,  a  practical 
and  efficient  channel  flow  routing  model  is  needed.  For  this  purpose,  the  Muskingum- 
Cunge  flow  routing  scheme  with  variable  parameters  (Koussis,  1978;  Laurenson,  1962; 
Ponce  and  Yevjevich,  1978)  has  been  modified  (Garbrecht  and  Brunner,  1991).  The 
hydrologic  approach  greatly  improves  computational  efficiency  and  speed,  and  reduces 
the  amount  and  detail  of  field  data  traditionally  needed  for  hydraulic  routing  (Weinmanr 
and  Laurenson,  1979).  Such  a  hydrologic  routing  scheme  is  a  practical  approach  to 
integrate  the  response  from  a  large  number  of  upstream  source  areas,  to  quantify  effects 
of  the  flow  integration  processes  on  watershed  runoff  characteristics,  and  to  investigate 
the  impact  of  spatially  variable  source-area  runoff  on  watershed  response. 

In  this  report,  the  hydrologic  Drainage  Network  Channel  Fir  /,  Routing  model 
(DNCFR)  is  presented,  followed  by  a  verification  for  channels  with  sample  and  compound 
sections,  and  a  third  order  drainage  network.  The  Muskingum-Cunge  flow  routing 
scheme  with  variable  parameters  is  used  as  the  initial  base  model.  It  is  adapted  for 
separate  flow  routing  in  the  main  and  overbank  charnel  portions,  and  it  includes  a 
variable  computational  time  increment.  The  parallel  main  and  overbank  channel  flow 
routing  simulates  the  flow  characteristics  in  each  channel  portion.  Differentiation  between 
main  and  overbank  channel  flow  is  often  desirable  because  sediment  mobilization, 
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transport  and  deposition,  transmission  losses  and  water  quality  parameters  vary  differently 
in  each  channel  portion  and  may  require  separate  treatment.  The  variable  computational 
time  increment  i"*  introduced  for  efficient  long-term  simulation  often  required  for  sediment 
and  water-quality  investigations. 

This  hydrologic  channel-flow  routing  scheme  is  applied  to  drainage  networks  by 
feeding  source-area  runoff  into  the  channels,  merging  the  channel  flows  at  network 
junctions,  and  routing  the  flows  through  the  channel  network.  As  for  most  hydrologic 
routing  schemes,  the  present  scheme  does  not  account  for  backwater  effects  and  does 
not  provide  detailed  hydraulic  flow  conditions  along  individual  channel  reaches  nor  does 
it  reproduce  localized  effects.  The  results  of  the  model  verification  show  that  the  DNCFR 
is  an  effective  tool  for  applications  to  large  complex  drainage  networks  and  for  continuous 
long-term  simulations.  To  operate  DNCFR  the  user  must  provide,  in  addition  to  the  usual 
channel  and  drainage  network  parameters,  surface  and  subsurface  inflows  into,  as  well 
as  losses  out  of  the  drainage  network. 

MODEL  DNCFR 

The  flow  routing  model  DNCFR  consists  of  four  components:  the  first  component 
quantifies  the  drainage  network  topology;  the  second  the  hydraulic  properties  of  the  cross 
sections;  the  third  routes  the  flow  in  individual  channel  reaches;  and  the  fourth  is  the  main 
driver  which  controls  the  execution  of  individual  model  components  and  coordinates  the 
routing  within  the  drainage  network.  Each  component  is  presented  separately. 

1  -  Drainage  Network  Topology  Component 

In  a  drainage  network,  it  is  generally  necessary  to  determine  the  sequence  in  which 
channel  flow  must  be  routed.  When  backwater  effects  are  negligible,  it  is  common 
practice  to  route  channel  flows  from  upstream  to  downstream.  Such  channel  flow  routing 
is  called  cascade  routing.  In  drainage  networks  there  are  many  upstream  channels  that 
simultaneously  contribute  to  the  watershed  outlet.  As  a  result,  there  are  many  possible 
sequences  in  which  channel  flows  can  be  routed.  The  determination  of  this  routing 
sequence  is  often  performed  manually.  This  is  a  tedious  and  error  prone  task  and  is  least 
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adaptable  to  subsequent  changes  in  drainage  network  resolution.  An  automated 
determination  of  the  routing  sequence  greatly  simplifies  the  engineer's  task,  insures 
correctness  and  consistency,  and  expedites  drainage  network  evaluation. 

One  such  programmable  algorithm  was  presented  by  Croley  (1980).  His  algorithm 
always  selects  the  right  most  source  node  as  starting  point,  and  it  determines  the 
sequence  of  channel  flow  routing  consistently  from  right  to  left  on  the  source  nodes, 
irrespective  of  drainage  network  configuration.  Under  certain  conditions  this  approach 
can  lead  to  large  computer  storage  needs,  as  subsequently  explained.  The  algorithm  of 
model  DNCFR  is  a  direct  solution  algorithm  (Garbrecht,  1988)  based  on  the  drainage 
network  definition  by  Croley  (1980).  It  determines  a  channel  flow  routing  sequence  that 
minimizes  computer  storage  needs.  Computer  storage  need  is  defined  as  the  number 
of  internal  arrays  required  for  cascade  routing.  An  array  is  necessary  for  temporary 
storage  of  runoff  results  from  one  channel  or  network  branch,  while  runoff  from  another 
is  being  evaluated.  For  example,  at  a  junction  node,  runoff  values  from  one  upstream 
inflow  must  be  stored,  while  the  other  upstream  inflow  is  being  evaluated.  This 
corresponds  to  one  storage  need.  Different  drainage  network  topologies  have  different 
computer  storage  needs.  In  the  following  the  algorithm  of  DNCFR  to  determine  the 
routing  sequence  is  briefly  presented. 

The  drainage  network  is  represented  as  an  arrangement  of  channels  and 
connecting  points  called  nodes.  The  type  of  node  is  defined  by  the  node  code.  A  node 
can  be  a  channel  source,  a  tributary  junction,  a  lateral  inflow  point,  or  any  other  special- 
purpose  point,  such  as  a  change  in  channel  cross  section  geometry  or  longitudinal  profile 
node.  A  list  of  node  codes  that  are  accepted  by  the  algorithm  is  given  in  Table  1. 
Between  two  nodes  channel  cross  section  geometry  and  longitudinal  slope  are  assumed 
constant.  It  is  also  assumed  that  junctions  of  more  than  two  channels  at  one  point  do 
not  occur.  However,  in  the  remote  chance  of  occurrence  the  situation  can  be  simulated 
by  adjacent  nodes  connected  by  a  short  channel.  Nodes  are  numbered  in  a  left  hand 
pattern  as  shown  in  Fig.  la  and  corresponding  node  codes  are  shown  in  Fig.  lb. 
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Table  1.  Node  code  definitions. 


Node  Code 

Definition 

1 

Source  node 

2 

Drainage  network  outlet  node 

3 

Channel  junction  node 

4 

Point  lateral  inflow  node 

5 

Water  withdrawal  node 

Change  in  channel  geometry  node 

Reservoir  inlet  node 

Reservoir  outlet  node 

Water  inflow  node  other  than  channel  junction  or  point  lateral  flow 

To  determine  the  optimal  channel  flow  routing  sequence  the  Strahler  channel 
orders  (Strahler,  1956)  at  each  node  is  needed.  This  is  accomplished  by  having  the 
algorithm  backtrack  from  the  source  nodes  downstream  and  increase  the  channel  order 
each  time  a  tributary  of  same  order  is  encountered.  When  a  tributary  is  of  large  Strahler 
order  then  the  latter  value  is  assumed  (Fig.  1c).  Once  Strahler  orders  are  assigned  to  all 
nodes,  upstream  and  lateral  inflows  to  each  node  must  be  identified.  The  algorithm 
identifies  upstream  inflows  into  a  junction  by  the  node  numbers  corresponding  to  the  two 
merging  channels,  and  lateral  inflow  by  the  node  number  it  flows  into. 

With  the  Strahler  channel  order  and  the  inflows  into  nodes  defined,  the  channel 
flow  routing  sequence  can  easily  be  determined.  The  flow  routing  begins  at  a  source 
node  that  leads  to  the  minimum  storage  needs.  This  source  node  is  one  that  directly 
contributes  to  the  network  order  (Garbrecht,  1988).  From  the  beginning  source  node,  the 
algorithm  backtracks  downstream  from  node  to  node,  and  from  network  subbranch  to  the 
next  larger  subbranch,  assigning  the  appropriate  routing  sequence  to  all  channels,  as 
shown  in  Fig.  Id.  Information  available  upon  completing  the  drainage  network  evaluation 
includes:  1)  channel  flow  execution  sequence,  2)  identification  of  upstream  and  lateral 
inflows,  3)  Strahler  channel  order  at  each  node,  4)  specification  of  channel  or  reservoir 
segments.  This  information  gives  a  complete  and  sufficient  description  of  the  drainage 
network  topology  to  fully  automate  the  management  of  the  channel  flow  routing  process. 


numbering  following  left 
er's  channel  orders;  d) 
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2  -  Cross  Section  Hydraulic  Properties  Component 

Hydraulic  properties  of  channel  cross-sections  (hereafter  referred  to  as  HPs)  are 
required  for  numerical  channel  flow  routing.  HPs  of  interest  are  cross-sectional  area,  top 
width  and  conveyance  factor.  They  are  a  function  of  stage,  and  therefore,  require 
repeated  evaluation  during  flow  routing  as  stage  varies  with  discharge.  This  calls  for  an 
efficient  scheme  to  quantify  the  HPs.  Model  DNCFR  uses  the  power  function  approach 
in  which  the  HPs  are  approximated  by  a  power  function  with  flow  depth  as  the 
independent  variable  (Li  et  al.,  1975;  Simons  et  al.,  1982;  Brown,  1982). 

HP  =  mD”  (1) 

where  HP  is  the  hydraulic  property,  D  is  flow  depth  and  m  and  p  are  coefficients  of  the 
power  function. 

The  coefficients  of  the  power  functions  are  computed  by  a  least  squared 
regression  through  the  logarithm  of  incremental  depth  and  HP  data  points.  This 
approximation  of  HPs  is  computationally  effective  and  generally  accurate  for  simple 
concave  sections.  In  the  case  of  compound  sections,  model  DNCFR  performs  the 
routing  separately  in  the  main  and  overbank  channel  portions,  as  previously  stated  and 
as  discussed  subsequently.  Therefore,  compound  sections  are  broken  into  two  simple 
sections,  and  two  separate  power  functions,  one  for  each  channel  portion,  are  developed 
as  illustrated  in  Fig.  2  for  the  wetted  perimeter.  It  is  assumed  that  the  power  function 
accurately  represents  the  rating  curve  for  simple  sections. 

3  -  Channel  Flow  Routing  Component 

The  channel  flow  routing  is  based  on  the  Muskingum-Cunge  routing  method  with 
variable  parameters,  with  further  adaptations  to  allow  for  variable  time  and  space 
increments,  and  routing  in  compound  sections.  Even  though  these  four  items  are  fully 
integrated,  they  are,  for  clarity  purposes,  presented  separately. 
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Figure  2.  Schematic  of  wetted  perimeter  representation  for  compound  cross  sections 

(not  to  scale). 
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-  Channel  Flow  Routing  Scheme 

The  channel  flow  routing  scheme  is  based  on  the  Muskingum-Cunge  routing 
method  with  variable  parameters.  The  Muskingum-Cunge  method,  and  refinements 
thereof,  have  been  amply  documented  in  previous  work  by  Cunge  (1969),  Kouscis  (1976, 
1980),  Ponce  (1983),  Ponce  and  Yevjevich  (1978),  Smith  (1980).  and  Weinmann  ana 
Laurenson  (1979).  The  method  is  a  kinematic  wave  routing  method.  The  kinematic  wave 
equation  is  transformed  into  a  diffusion  equation  by  numerical  attenuation  of  the 
imperfectly  centered  finite-difference  scheme  (Smith,  1980).  The  method  therefore 
accounts  for  hydrograph  convection  and  diffusion,  i.e.  for  downstream  movement  and 
peak  attenuation  of  the  hydrograph.  Diffusion  is  introduced  through  two  weighting 
coefficients  which  are  determined  from  physical  channel  properties  and  flow 
characteristics  (Cunge,  1969).  When  these  coefficients  are  varied  as  a  function  of  flow 
the  method  becomes  a  non-linear  coefficient  method  (Koussis,  1976;  Laurenson,  1962). 
The  Muskingum-Cunge  routing  method  with  variable  parameters  accounts  for  most  of  the 
flood  wave  phenomena  when  practical  applications  are  considered  (Ponce  and  Theurer, 
1982;  Weinmanr  and  Laurenson,  1979).  The  advantages  ot  this  method  over  other 
hydrologic  techr.'':;:.es  such  as  normal  depth.  Modified  Puls,  or  simple  Muskingum  method 
are:  (1)  the  scheme  is  stable  with  properly  selected  coefficients  (Smith,  1980;  Ponce, 
1981;  Ponce  and  Theurer,  1982);  (2)  it  produces  consistent  results  in  that  the  results  are 
reproducible  with  varying  grid  resolution  (Jones,  1983;  Koussis,  1983;  Ponce  and  Theurer, 
1982,  1983a  and  1983b);  (3)  it  is  comparable  to  the  diffusion  wave  routing  (Cunge,  1969; 
Miller  and  Cunge,  1975);  (4)  the  coefficients  of  the  method  are  physically  based  (Cunge, 
1969);  (5)  the  method  has  been  shown  to  compare  well  against  the  full  unsteady  flow 
equation  over  a  wide  range  of  flow  situations  (Ponce,  1981;  Younkin  and  Merkel,  1988a 
and  1988b);  and  (6)  the  solution  is  largely  independent  of  time  and  space  intervals  when 
these  are  selected  within  the  spatial  and  temporal  resolution  criteria  (Ponce,  1981 ;  Ponce 
and  Theurer,  1982).  The  essential  steps  of  this  method  are  briefly  summarized  in  the 
following.  Detailed  formulations  and  discussions  can  be  found  in  the  above  cited 
literature. 
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The  Muskingum-Cunge  routing  schenrie  uses  a  storage  relation  to  relate  inflow  and 
outflow  in  a  channel  reach.  The  storage  relation  is  given  by: 

S  =  K[XI  +  (1  -  X)  O]  (2) 

where  K  is  a  storage  coefficient,  X  is  a  weighting  factor,  I  is  the  inflow  rate  to  the  reach, 
and  O  is  the  outflow  rate  from  the  reach.  The  finite  difference  formulation  of  Eq.  2  results 
in  the  Muskingum  Equation  (Cunge,  1969;  Weinmann  and  Laurenson,  1979): 


with 


n  +  l  n  n  +  1  n 

O  .  =  C  Q  .  +  C  Q  .  +CO. 

+  l  j  +  1 


(3) 


At 


+  2X 


= 


(4) 


At 


-  2X 


= 


(5) 


^3  = 


At 

2  (1  -  X)  -  — 
K 


(6) 


At 

F  =  -  +  2  (1  -  X)  (7) 

K 

where  n  is  time  superscript,  j  is  space  subscript,  Q  is  discharge,  and  At  is  the  routing  time 
increment  of  the  finite  difference  cell.  In  the  original  Muskingum  equation,  the  value  of  the 
storage  coefficient  K  and  the  weighting  factor  X  are  determined  by  trial  and  error  or  by 
calibration  with  observed  hydrographs  (Miller  and  Cunge,  1975).  In  the 
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Muskingum-Cunge  approach  coefficients  K  and  X  are  expressed  in  terms  of  flow,  channel 
and  finite  difference  cell  parameters  (Cunge,  1969;  Koussis,  1978;  Ponce  and  Theurer, 
1982;  Weinmann  and  Lawrence,  1979)  as: 

Ax 

K  =  -  (8) 

c 


where  ax  is  the  space  increment  of  the  finite  difference  cell,  c  is  a  representative 
floodwave  celerity,  q  is  a  representative  unit  width  discharge,  and  is  the  channel  bed 
slope.  With  Eq.  8  and  9,  the  need  of  observed  hydrographs  to  calibrate  the  coefficients 
K  and  X  is  eliminated.  Cunge  (1969)  also  demonstrated  that  the  Muskingum-Cunge 
scheme,  given  by  Eqs.  3  through  9,  is  equivalent  to  a  convection-diffusion  wave  model, 
i.e.  accounting  for  downstream  movement  and  peak  attenuation  of  the  hydrograph. 

Discharge  and  flood  wave  celerity  are  generally  different  at  various  points  along  a 
flood  wave.  To  account  for  some  of  this  observed  nonlinearity,  Koussis  (1976), 
Weinmann  and  Laurenson  (1980),  and  Ponce  and  Yevjevich  (1978)  presented  the  concept 
of  variable  coefficients.  They  redefined  coefficients  K  and  X  for  every  computational  cell 
as  a  function  of  updated  values  of  unit  width  discharge  and  wave  celerity.  The  unit  width 
discharge  and  wave  celerity  at  a  grid  point  (j,  n)  are  defined  as: 


(10) 


(11) 


where  Q  is  total  discharge,  A  is  flow  area,  B  is  top  width,  and  c  is  the  floodwave  celerity. 
The  celerity  is  derived  from  the  equation  of  continuity  following  the  Kleitz-Seddon  principle 
(Chow,  1959).  The  relation  between  discharge  and  flow  area  (Eq.  10)  is  based  on 
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Manning’s  uniform  flow  equation  with  energy  slope  equal  to  bed  slope.  It  is  therefore  a 
kinematic  wave  celerity.  The  average  flood  wave  celerity  for  a  computational  cell  is  given 
as  the  average  value  of  the  celerity  at  the  four  nodes  of  the  cell. 

For  a  computational  cell,  the  unknown  unit  width  discharge  and  wave  celerity  are 
evaluated  by  a  four-point  iterative  approximation  (Ponce  and  Yevjevich,  1978).  To  begin 
the  iteration  an  initial  estimate  of  the  discharge  for  the  unknown  grid  point  (j+1.  n  +  1)  is 
obtained  using  a  linear  projection  of  the  known  discharge  at  points  (j,  n),  (j+1,  n)  and 
(j,  n+1).  Thereafter,  a  four-point  iteration  is  used  to  solve  for  the  discharge  at  the 
unknown  point.  The  relation  between  discharge,  flow  area,  top  width,  and  flow  depth  is 
defined  by  power  functions  which  are  derived  using  cross  section  shape  and  Manning's 
uniform  flow  equation.  These  power  functions  represent  simple  rating  curves. 

The  Muskingum-Cunge  method  with  variable  parameters  was  found  to  be  accurate 
for  a  wide  range  of  simple  channels  and  flow  conditions  (Younkin  and  Merkel  (1988a  and 
1988b)).  Younkin  and  Merkel  (1988b)  performed  340  routing  tests  and  compared  the 
results  to  those  from  a  full  dynamic  model  used  as  a  benchmark.  They  found  that  peak 
discharge,  peak  area,  times  to  peak,  and  correlation  of  hydrograph  shapes  satisfied  over 
80  per  cent  of  the  Soil  Conservation  Service  (SCS)  field  conditions  covered  by  their  study. 
The  accuracy  criteria  used  in  their  study  were:  (1)  less  than  one  percent  difference  for 
peak  discharge  and  area;  (2)  one  or  less  time  step  difference  between  time  of  occurrence 
of  discharge  and  area  peaks;  and,  (3)  greater  than  ninety -five  percent  shape  correlation 
for  discharge  and  area  hydrographs.  However,  the  flow  routing  scheme  does  not 
account  for  backwater  effects  and  it  diverges  from  the  full  unsteady  flow  solution  for  very 
rapidly  rising  hydrographs  in  flat  channels  with  slope  less  than  0.0002  (Brunner,  1989). 

-  Variable  Computational  Time  Increment 

Variable  computational  time  increments  are  introduced  to  increase  numerical 
efficiency  of  the  routing  scheme.  Large  time  increments  are  used  during  inter-storm 
periods  when  relatively  constant  discharges  prevail.  Shorter  time  increments  apply  when 
discharge  varies  rapidly  during  rainfall-runoff  events.  The  change  in  time  increment  is 
gradual  to  assure  smooth  transitions  and  to  prevent  a  hydrograph  from  moving  from  a 
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region  with  fine  computational  cells  to  one  with  coarse  ones. 

Variable  time  increments  are  compatible  with  the  finite  difference  formulation  of  the 
Muskingum-Cunge  routing  scheme  (Eqs.  3  through  7)  because  the  latter  is  an  explicit 
scheme.  Flow  calculations  dependent  only  on  the  current  space  and  time  increment,  and 
they  are  independent  of  any  other  computational  cell.  As  a  result  there  are  no 
requirements  to  keep  ax  and  At  constant  throughout  the  computational  domain,  and  they 
may  vary  within  limits  established  by  the  accuracy  criteria  set  forth  by  Ponce  and  Theurer 
(1982). 

The  size  of  the  time  increment  is  determined  as  a  function  of  rate  of  change  in 
upstream  inflow  into  the  channel  reach.  The  inflow  time  series  is  scanned  ahead.  If  a 
change  in  discharge  above  a  given  threshold  value  is  sensed,  the  current  time  increment 
is  reduced;  if  no  change  in  discharge  is  found,  the  size  of  the  next  time  increment  is 
increased.  Upper  and  lower  bounds  for  the  time  increment  size  are  one  day  and  five 
minutes,  respectively.  These  boundaries  were  found  to  work  well  for  long-term  simulation 
in  drainage  networks.  In  addition  to  the  smooth  transition  between  fine  and  coarse  time 
increments,  the  early  reduction  of  the  time  increment  size  as  an  upcoming  perturbation 
is  sensed  assures  an  adequate  temporal  resolution  for  hydrograph  routing  that  generally 
satisfies  the  accuracy  criteria  of  a  minimum  of  5  time  increments  on  the  rising  portion  of 
a  hydrograph  (Ponce  and  Theurer,  1982). 

-  Computational  Space  Increment 

Computational  space  increments,  A  x,  are  subreaches  that  define  the  computational 
cell  size  at  which  the  numerical  flow  routing  is  performed.  A  computational  space 
increment  may  be  equal  to  the  entire  routing  reach  length  or  to  a  fraction  of  that  length. 
It  is  initially  selected  as  the  entire  reach  length.  If  the  size  of  this  space  increment  does 
not  meet  the  accuracy  criteria  for  flow  routing  given  by  Ponce  and  Theurer  (1982),  it  is 
reevaluated  by  subdividing  the  length  of  the  routing  reach  into  even  subreaches  that 
produce  A  x’s  that  satisfy  the  accuracy  criteria.  Ponce  and  Theurer's  accuracy  criteria  is 
given  by: 
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with 


AJf  ^d) 


X  =  2 


(12) 


(13a) 


Ax  =  cAt  (13b) 


where  q  and  c  refer  to  a  reference  discharge  and  celerity,  respectively,  X  is  an  accuracy 
parameter  and  At  is  the  minimum  time  increment.  The  minimum  time  increment  is  used 
because  it  is  the  one  applicable  during  routing  of  a  hydrograph.  The  reference  discharge 
is  generally  two  thirds  of  the  peak  flow  above  base  flow,  and  the  reference  celerity  is  the 
celerity  corresponding  to  the  reference  discharge. 

The  upper  limit  of  the  space  increment,  as  given  by  Eq.  12,  becomes  quite  large 
for  very  flat  channels  and  high  discharge  values.  In  long  channel  reaches  where  such 
large  space  increments  can  be  implemented,  the  flow  routing  may  produce  inaccurate 
hydrographs.  First,  the  time  separation  between  inflow  and  outflow  hydrographs  can 
become  large  resulting  in  the  computed  outflow  hydrograph  to  end  up  in  a  region  of 
coarse  time  increments.  In  this  case  the  upper  limit  of  the  space  increment  depends  on 
the  duration  and  celerity  of  the  hydrograph.  Short  duration  and  fast  moving  hydrographs 
require  shorter  space  increments  than  long  duration  and  slow  moving  hydrographs.  As 
a  rule  of  thumb  the  average  hydrograph  travel  time  in  a  space  increment  should  not 
exceed  about  one  fifth  of  the  duration  of  the  inflow  hydrograph. 

In  the  second  case,  during  overbank  flow  conditions,  long  space  increments  may 
result  in  the  hydrograph  in  the  main  channel  to  significantly  outpace  the  hydrograph  on 
the  overbank  portion  of  the  cross  section.  This  outpacing  is  a  natural  phenomena  that 
changes  hydrograph  shape.  Flow  from  main  channel  hydrograph  spills  onto  the  flood 
plains,  peak  runoff  rate  decreases,  and  the  recession  limb  is  stretched  out  as  a  result  of 
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return  flow  from  the  overbanks.  In  the  present  routing  scheme,  the  uncoupling  of  the 
main  channel  and  overbank  flow  routing  in  very  long  space  increments  may  produce  an 
insufficient  flow  mixing  between  the  two  channel  portions  and,  as  a  result,  the  effect  of 
flood  plain  storage  on  hydrograph  shape  is  not  accurately  simulated.  An  upper  limit  for 
the  space  increment  of  about  one  twentieth  of  the  wave  length  was  found  to  provide,  in 
most  cases,  an  adequate  flow  mixing.  Under  real  world  conditions  tributary  junctions  in 
drainage  networks  and  changes  in  cross  section  and  flood  plain  characteristics  generally 
provide  short  enough  channel  reach  length  that  do  not  require  limitations  on  the  space 
increment. 

-  Routing  in  Compound  Cross  Sections 

Main  and  overbank  channel  portions  are  separated  and  modeled  as  two 
independent  channels.  Right  and  left  overbanks  are  combined  into  a  single  overbank 
channel.  At  the  upstream  end  of  a  space  increment  total  inflow  discharge  is  divided  into 
a  main  channel  and  an  overbank  flow  component.  Each  is  then  routed  independently 
using  the  previously  described  routing  scheme.  Momentum  exchange  at  the  flow 
interface  between  the  two  channel  portions  is  neglected  and  the  hydraulic  flow 
characteristics  are  determined  for  each  channel  portion  separately.  At  the  downstream 
end  of  the  space  increment  both  flow  components  are  summed  to  yield  the  total  outflow 
discharge.  The  flow  exchange  between  main  channel  and  overbank  channel  during 
routing  within  a  space  increment  is  neglected. 

Flow  redistribution  between  main  and  overbank  channels  is  based  on  the 
assumption  of  a  constant  energy  head  perpendicular  to  the  flow  direction  and  on  a 
negligible  momentum  exchange  at  the  flow  interface  between  the  two  channel  portions. 
As  the  stage  in  the  main  channel  exceeds  overbank  elevation,  the  discharge  in  each 
channel  portion  is  determined  by  matching  the  energy  head  of  the  flow.  The  energy  head 
is  computed  using  Bernoulli's  conservation  of  energy  equation  and  mean  flow  values  for 
each  channel  portion: 
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Z  d  +  —  =  Z  d  +  —  (14) 

2g  °  °  2g 

0,  =  (15) 

with  V  =  Q/A;  A  =  fct  (Q,  XSS,  N);  and  d  =  fct  (Q,  XSS,  N);  where  Z  is  elevation  above 
a  reference  datum,  v  is  flow  velocity,  XSS  is  channel  section  shape,  N  is  channel 
roughness  and  d  is  flow  depth.  Subscript  m  stands  for  main  channel,  o  for  overbank 
channel,  and  T  for  total.  Flow  area.  A,  and  flow  depth,  d,  are  determined  as  a  function 
of  cross  section  shape  and  Manning's  uniform  flow  equation.  The  described  separation 
into  main  and  overbank  flow  components  introduces  a  lateral  variation  in  flow 
characteristics  which  make  the  routing  scheme  a  quasi  two  dimensional  approach.  From 
the  point  of  view  of  flow  routing  the  uncoupling  of  main  and  overbank  flow  results  in  a 
flood  wave  propagation  that  is  primarily  controlled  by  the  faster  moving  flow  in  the  main 
channel,  and  a  wave  attenuation  that  is  primarily  controlled  by  the  storage  of  the 
overbank  channel. 

4  -  Coordination  Component  of  Drainage  Network  Routing 

This  model  component  uses  the  network  topology  data  determined  in  model 
component  1  to  coordinate  the  routing  in  the  drainage  network.  It  defines  and  feeds  the 
source  area  runoff  into  the  channels,  merges  appropriate  channel  flows  at  network 
junctions,  and  executes  the  flow  routing  in  proper  sequence  for  all  channel  reaches. 
Because  this  model  component  performs  simple  bookkeeping  tasks,  no  further 
explanations  are  given. 

VERIFICATION  APPROACH 

The  original  Muskingum-Cunge  channel  flow  routing  with  variable  parameters  was 
found  to  be  accurate  for  a  wide  range  of  channel  geometries  and  flow  conditions 
(Koussis,  1978;  Ponce,  1981;  Ponce  and  Theurer,  1982;  Younkin  and  Merkel,  1988a  and 
1988b).  Brunner  (1989)  indicated  that  the  method  diverges  from  the  full  unsteady  flow 
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solution  for  very  rapidly  rising  hydrographs  in  flat  channel  with  slopes  less  than  0.0002. 
In  the  context  of  this  report,  DNCFR  is  tested  for  channels  with  simple  and  compound 
sections  and  for  complex  drainage  networks. 

Verification  is  accomplished  by  comparing  flow  routing  results  with  corresponding 
results  from  models  solving  the  full  unsteady  flow  equations.  This  benchmark  verification 
approach  is  preferred  over  actual  field  data,  because  field  data  often  includes  processes 
unrelated  to  flow  routing  such  as  infiltration  on  overbank  flood  plains,  or  variable  flow 
resistance  due  to  the  submergence  of  vegetation.  These  effects  are  generally  hard  to 
measure,  and  calibration  to  site-specific  conditions  is  often  required.  The  latter  makes 
any  comparison  between  model  results  and  field  conditions  highly  subjective.  Benchmark 
verification  assures  well  defined  and  identical  boundary  conditions.  It  provides  an 
objective  comparison  to  state-of-the-art  one-dimensional  hydraulic  modeling  capabilities. 

The  models  selected  as  benchmark  are  DAMBRK  (1988  version)  of  the  National 
Weather  Service  (Fread,  1984)  and  UNET  (Version  1.1)  by  Barkau  (1990).  The  DAMBRK 
model  is  used  for  verification  of  flow  routing  in  single  channels  with  simple  and  compound 
sections.  The  UNET  model  is  used  to  verify  the  flow  routing  in  drainage  networks.  These 
models  were  selected  because  of  the  solution  to  the  full  dynamic  flow  equations  and  no 
model  comparison  is  intended. 

Verification  criteria  are  peak  discharge,  time  to  peak,  runoff  volume  and 
hydrograph  shape.  Discrepancies  from  the  benchmark  are  quantified  in  percent  deviation 
for  the  first  three  parameters.  Hydrograph  shape  is  verified  visually  and  quantified  by  the 
correlation  coefficient  between  computed  and  benchmark  hydrograph  values.  Finally,  the 
numerical  performance  of  DNCFR  versus  the  hydraulic  benchmark  models  is  defined  in 
terms  of  a  reduction  in  overall  execution  time. 

VERIFICATION  TEST  CASES 

Eighteen  verification  cases  are  presented.  Of  these,  ten  are  single  channels  with 
simple  cross  sections,  six  are  single  channels  with  compound  sections,  and  two  are 
drainage  networks  with  a  mix  of  channels  with  simple  and  compound  sections.  Channel 
geometry,  resistance  to  flow,  and  inflow  hydrograph  characteristics  for  the  single  channel 


19 


tests  are  given  in  Table  2.  A  schematic  of  the  trapezoidal  channel  cross  sections  for  test 
cases  2.5  and  2.6  are  shown  in  Fig.  3.  The  inflow  hydrographs  are  generated  using  a 
gamma  function  (Ponce  and  Theurer,  1982).  The  inflow  hydrograph  of  test  1.4  has  a 
dimensionless  wave  period,  t  ,  greater  than  171  (Table  2),  and,  according  to  Ponce  at  al. 
(1978),  it  is  classified  as  a  kinematic  wave.  Of  the  other  inflow  hydrographs  seven  are 
diffusion  waves  with  a  t  /F,  factor  greater  than  30,  and  eight  are  dynamic  waves  with  a 
T  /F,  factor  under  30.  The  predominance  of  diffusion  and  dynamic  inflow  hydrographs 
makes  the  selected  hydrographs  relevant  for  the  testing  of  DNCFR. 


TEST  CASE  2.5 
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Figure  3.  Schematic  of  cross  section  for  test  cases  2.5  and  2.6  (distorted  vertical  scale). 

The  hypothetical  drainage  network  for  the  two  network  tests  are  shown  in  Fig.  4. 
The  watershed  is  approximately  17  km  long  and  9  km  wide.  Channel  slopes  vary  from 
0.0016  to  0.004;  first  and  second  order  channels  have  no  overbank  flood  plains;  third 
order  channels  have  significant  flood  plains.  The  inflow  hydrographs  are  mostly  of  the 
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Figure  4.  Hypothetical  drainage  network  for  tests  3.1  and  3.2. 

diffusion  type,  and  are  entered  at  source  nodes  and  selected  junction  nodes  as  indicated 
by  hollow  circles  in  Fig.  4.  Peak  and  time  to  peak  of  the  inflow  hydrograph  are  chosen 
arbitrarily;  they  range  from  3.0  to  7.5  cms,  and  55  to  80  min,  respectively,  for  test  3. 1 ;  and 
from  3.0  to  6.5  cms  and  145  to  170  min,  respectively,  for  test  3.2.  Hydrograph  shapes 
are  generated  using  a  gamma  function.  A  storm  movement  at  about  20  km/hr  is 
assumed  from  the  outlet  to  the  top  of  the  basin. 

RESULTS  OF  VERIFICATION 

The  results  of  the  verification  are  shown  in  Table  3.  For  channels  with  simple 
sections,  hydrograph  peak  and  time  to  peak  are,  on  the  average,  less  than  3%  off  the 


Table  3:  Verification  results  of  benchmark  and  DNCFR,  and  differences  between  the  two. 
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Duration  of  simulation:  24  hrs;  minimum  time  increments;  5  min.;  includes  I/O. 
NO:  No  data;  NA:  not  applicable. 


benchmark  values;  for  compound  sections  they  are,  on  the  average,  less  than  4%  off  the 
benchmark  values.  In  general,  tests  having  hydrographs  of  the  dynamic  type  (tests  1.5, 
1.6,  1.7,  1.8,  2.1,  2.3,  and  2.4)  display  a  larger  discrepancy  than  those  having 
hydrographs  of  the  diffusion  type  (tests  1.1,  1.2,  1.3,  2  2,  2.5,  and  2.6).  This  is  to  be 
expected  because  DNCFR  is  equivalent  to  a  diffusion  routing  method  (Cunge,  1969),  and 
it  does  not  account  for  dynamic  effects.  Considering  the  highly  dynamic  character  of 
some  of  the  hydrographs,  the  results  of  DNCFR  are  good  for  a  hydrologic  routing 
method.  For  example,  in  channels  with  simple  sections  and  dynamic  hydrographs  (tests 
1.5,  1.6,  1.7,  and  1.8),  about  91%  of  the  attenuation  is  reproduced  by  DNCFR.  For 
channels  with  compound  sections,  about  94%  of  the  total  attenuation  (due  to  diffusion 
and  storage)  is  reproduced  by  DNCFR. 

Figure  5  show  a  plot  of  discharge  and  time  to  peak  from  DNCFR  versus 
corresponding  benchmark  values.  The  data  represents  all  tests  of  Table  1 ,  including  the 
two  network  applications.  Neither  discharge  nor  the  time  to  peak  show  significant 
deviations  from  the  line  of  perfect  agreement.  In  the  following  peak  discharge,  time  to 
peak  and  hydrograph  shape  are  discussed  in  more  detail  for  selected  tests. 

With  respect  to  the  drainage  netv/ork  tests,  the  complex  hydrographs  are  the  result 
of  the  drainage  network  configuration  and  movement  of  the  storm  up  the  watershed.  The 
runoff  from  the  lower  right  network  branch  arrives  first  at  the  outlet  follow  ^d  by  the  main 
peak  from  the  upper  portion  of  the  watershed.  The  higher  peak  values  by  DNCFR  are 
primarily  the  result  of  limitations  regarding  backwater  and  reverse  flow  effects.  Indeed, 
the  hydraulic  simulations  include  reverse  flow  up  tributary  branches.  Reverse  flow  is  the 
result  of  high  stages  in  the  main  channel  while  low  stages  prevail  on  the  tributaries.  The 
net  effect  of  this  reverse  flow  is  additional  storage  and  attenuation  of  the  peak.  Times  to 
peak  and  hydrograph  shapes  are  well  reproduced  for  both  network  applications.  The 
hydrographs  leaving  the  drainage  network  are  shown  in  Fig.  6.  Peak,  timing  and  shapes 
are  well  reproduced.  The  correlation  coefficient  for  hydrograph  shape  is  0.98. 

For  simple  and  compound  sections,  the  average  correlation  coefficient  for 
hydrograph  shape  is  0.99  and  0.97,  respectively.  Test  1 .6  is  a  dynamic  wave  with  a  35% 
attenuation  of  the  peak  above  base  flow.  Ir^  this  case,  like  in  all  other  dynamic  cases,  the 
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Figure  5.  Benchmark  vs.  computed  discharge  and  time  to  peak;  all  verification  test 
results. 
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Figure  6.  Hydrographs  for  tests  3.1  and  3.2. 
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hydrograph  computed  by  DNCFR  displays  a  larger  attenuation  than  the  benchmark.  As 
a  result  of  the  lower  flow,  the  time  to  peak  lags  slightly  behind  the  benchmark  value.  Test 
1 .9  is  a  diffusion  wave  with  11%  attenuation  of  the  peak  above  base  flow.  It  is  typical  of 
most  diffusion  wave  cases  with  little  discrepancy  in  peak,  time  to  peak  and  hydrograph 
shape. 

Outflow  hydrographs  for  channels  with  compound  section  are  given  by  tests  2.1, 
2.2,  2.5  and  2.6  which  are  shown  in  Fig.  7.  Test  2.1  shows  a  hydrograph  that  is 
attenuated  entirely  below  overbank  flood  plain  elevation.  This  example  also  shows  the 
effect  of  the  slower  moving  overbank  flow  by  producing  a  longer  hydrograph  recession 
limb.  In  the  other  three  test  cases,  overbank  flow  is  active  along  the  entire  channel 
length.  For  these  cases  the  location  of  beginning  of  overbank  flow  is  clearly  defined  by 
the  breaks  in  the  rising  limb  of  the  hydrographs.  Overall,  the  distortions  of  hydrograph 
shape  due  to  overbank  flood  plains  are  consistent  with  the  benchmark  shapes. 

Finally,  the  execution  time  of  DNCFR  is,  on  the  average,  92%  shorter  than  for  the 
benchmark  hydraulic  evaluation  (Table  3).  The  comparison  is  made  on  an  IBM 
compatible  PC,  having  an  80836  -  20MHZ  micro-processor  and  math  co-processor,  and 
using  the  Microsoft  FORTRAN  compiler  Version  5.0  with  optimization.  Program  I/O  and 
support  computations  are  included  in  the  execution  time.  Using  the  average  execution 
time  over  all  tests,  the  reduction  is  from  10  minutes  to  32  seconds,  which  is  a  factor  of 
20,  or  about  1 .2  orders  of  magnitude.  It  is  believed  that  additional  reductions  in  execution 
time  can  be  achieved  for  long-term  simulations  because  significantly  larger  time 
increments  are  used  by  the  model  for  nearly  constant  discharge  values  that  generally 
prevail  between  storm  events. 
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SUMMARY  AND  CONCLUSIONS 

The  Muskingum-Cunge  channel  flow  routing  scheme  with  variable  parameters  is 
modified  to  account  for  compound  sections,  to  include  a  variable  time  step,  and  to 
determine  internally  the  computational  reach  increment.  The  resulting  model,  DNCFR,  is 
verified  for  hypothetical  channel  and  flow  conditions.  Routing  results  are  compared  with 
those  of  hydraulic  models  solving  the  full  unsteady  flow  equations.  Ten  channels  with 
simple  sections,  six  with  compound  sections  and  two  drainage  network  applications  are 
selected  for  verification.  For  ail  tested  cases,  DNCFR  reproduces  the  peak,  time  to  peak 
and  shape  of  the  benchmark  hydrograph  with  reasonable  accuracy.  Slight  discrepancies 
(less  than  10%)  in  the  drainage  network  application  are  due  to  the  limitations  of  hydrologic 
models  with  respect  to  backwater  and  reverse  flow  effects.  The  size  of  the  discrepancy 
is  well  within  the  usual  error  of  drainage  network  parameterization  and  lateral  channel 
inflow  determination.  Program  efficiency,  as  measured  by  the  reduction  in  execution  time, 
is,  on  the  average,  one  order  of  magnitude  faster  than  the  benchmark  hydraulic  routing. 
The  results  of  the  verification  indicate  DNCFR  to  be  an  effective  tool  for  hydrologic  routing 
applications  to  large  complex  drainage  networks  and  for  continuous  long-term 
simulations. 
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